package twoD.IsoPfem.Plasticity.J2PlaneStress;

import inf.jlinalg.BLAM;
import inf.jlinalg.IMatrix;
import inf.jlinalg.IVector;
import oneD.fem.algorithm.ITensor;
import oneD.fem.algorithm.Tensor2;
import oneD.fem.algorithm.Tensor4;
import twoD.IsoPfem.Plasticity.IHardeningLaw;
import twoD.IsoPfem.Plasticity.IYieldCriterion;

public class vonMisesYieldCriterion implements IYieldCriterion {

	private double alpha_;
	private Tensor2 beta_;
	private IHardeningLaw isotropicHardeningLaw_;
	private IHardeningLaw kinematicHardeningLaw_;

	public vonMisesYieldCriterion(IHardeningLaw iLaw, IHardeningLaw kLaw) {
		isotropicHardeningLaw_ = iLaw;
		kinematicHardeningLaw_ = kLaw;
	}

	public void setInternalVariables(double alpha, Tensor2 beta) {
		alpha_ = alpha;
		beta_ = beta;
	}

	@SuppressWarnings("unused")
	private double computeJ2(IVector sigma) {
		Tensor2 Sigma = new Tensor2(sigma);
		Tensor2 S = Sigma.getDeviatoricTensor();
		ITensor S2 = S.multiply(S);
		return (Math.pow(S.trace(), 2) - S2.trace()) / 2;
	}

	@Override
	public double f(IVector sigma) {
		Tensor2 Sigma = new Tensor2(sigma);
		Tensor2 S = Sigma.getDeviatoricTensor();
		IMatrix dev_sigma = S.getA();
		BLAM.add(-1.0, beta_.getA(), dev_sigma);
		double tmp = BLAM.norm(BLAM.NORM_FROBENIUS, dev_sigma);
		tmp -= Math.sqrt(2.0 / 3) * isotropicHardeningLaw_.valueAt(alpha_);
		return 0;
	}

	public double f(ITensor xi, double alpha) {
		double tmp = Math.sqrt(xi.dot(xi));
		tmp -= Math.sqrt(2.0 / 3) * kinematicHardeningLaw_.valueAt(alpha);
		return tmp;
	}

	public Tensor2 dfdsigma(IVector sigma) {
		Tensor2 Sigma = new Tensor2(sigma);
		Tensor2 S = Sigma.getDeviatoricTensor();
		IMatrix dev_sigma = S.getA();
		BLAM.add(-1.0, beta_.getA(), dev_sigma);
		double tmp = BLAM.norm(BLAM.NORM_FROBENIUS, dev_sigma);
		BLAM.scale(1 / tmp, dev_sigma);
		return new Tensor2(dev_sigma);
	}

	public double dfdalpha(IVector sigma) {
		return Math.sqrt(2.0 / 3);
	}

	public Tensor2 dfdbeta(IVector sigma) {
		Tensor2 N = dfdsigma(sigma);
		BLAM.scale(2.0 / 3 * kinematicHardeningLaw_.valueAt(alpha_), N.getA());
		return N;
	}

	public IHardeningLaw getIsotropicHardeningLaw() {
		return isotropicHardeningLaw_;
	}

	public IHardeningLaw getKinematicHardeningLaw() {
		return kinematicHardeningLaw_;
	}

	@Override
	public Tensor4 computeConstitutiveTangentOperator() {
		// TODO Auto-generated method stub
		return null;
	}
}
